home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-04 / extenpr.zip / EXTEND.FOR < prev    next >
Text File  |  1988-06-04  |  16KB  |  709 lines
  1. C
  2. C    This is a collection of subroutines that perform extended precision
  3. C    floating point arithmetic.
  4. C
  5. C    By  Mark P. Esplin
  6. C        16 Shawsheen Rd.
  7. C        Billerica, Ma 01821    
  8. C
  9. C    Whole programs do not have to be written using extended
  10. C    precision.  Only the part that needs extra precision need use it,
  11. C    then convert back to regular REAL*8 numbers.
  12. C    The routines for addition and multiplication use the same algorithms 
  13. C    that are taught in grade school but using base 10000 numbers.
  14. C    This is not the most  efficient way to do it, but
  15. C    it makes it possible to use standard FORTRAN77 (almost).  
  16. C    These routines should work with any computer and compiler that supports
  17. C    INTEGER*2 and INTEGER*4 arithmetic.
  18. C    Each number is stored in an INTEGER*2 array of length N.
  19. C    The first word, IA(1), is the exponent in base 10000 digits.
  20. C    The least significant digit is in IA(2) and the most significant
  21. C    digit is in IA(N).  If IA(N) > 5000 the number is negative.
  22. C    When IA(N) = 5000 the number is assumed invalid.
  23. C    Base 10000 digits have worse roundoff characteristics than binary.
  24. C    Since roundoff causes the last base 10000 digit to be inaccurate, 
  25. C    that really means the last 4 decimal digits are inaccurate.  Be
  26. C    sure to use more digits than is necessary to allow for roundoff.
  27. C    N must be at least 6 or some routines won't work correctly.
  28. C
  29. C    Subroutines with their FORTRAN equivalent in this package are:
  30. C
  31. C    ADD(A,B,C,N)      C=A+B
  32. C    NEG(A,N)          A=-A
  33. C    MUL(A,B,C,N)      C=A*B
  34. C    INV(A,B,N,W)      B=1/A      W is work area of size 3*N
  35. C    EXSQRT(A,B,N,W)   B=SQRT(A)  W is work area of size 4*N
  36. C    EXHALF(A,N)       A=1/2      much faster than using INV
  37. C    PRTEXT(A,N)       WRITE(*,*) A
  38. C    WRTEXT(IU,A,N)    WRITE(IU,*) A
  39. C    CVI2EX(I2,A,N)    A=I2       I2 is INTEGER*2
  40. C    CVR8EX(R,A,N)     A=R        R is REAL*8
  41. C    CVEXR8(A,N,R)     R=A        R is REAL*8
  42. C    COPY(A,B,N)       B=A
  43. C    IEXCMP(A,B,N,W)    Compares two positive numbers A and B returns 1, 0, -1
  44. C    CHGPRE(A,NA,B,NB)  Changes the precision of A form NA to NB
  45. C
  46. C    There are other internal routines, but they would probably be of limited
  47. C    value for most users.
  48. C
  49. C
  50.     SUBROUTINE INCEX(IA,N)
  51.     INTEGER*2 IA(*),N
  52. C    This routine increments the mantissa of extended number IA.
  53. C    This routine is used by other routines to perform rounding.
  54.     INTEGER*2 IT,ICAR,I
  55.     ICAR=1
  56.     DO 10 I=2,N
  57.       IT=IA(I)+ICAR
  58.       ICAR=IT/10000
  59.       IA(I)=IT-ICAR*10000
  60. C    In almost all cases this routine will end after a few passes.
  61.     IF(ICAR.EQ.0) GOTO 20
  62. 10    CONTINUE
  63. C    Shift back if necessary
  64. C    The only time that a shift is necessary is when IA(N) = 5000 and
  65. C    the rest of the mantissa is zero.
  66. 20    CONTINUE
  67.     IF(IA(N).EQ.5000) THEN
  68.       IA(1)=IA(1)+1
  69.       DO 30 I=2,N-1
  70.         IA(I)=IA(I+1)
  71. 30      CONTINUE
  72.       IA(N)=0
  73. C      Don't need to round again since if a shift was necessary
  74. C      the least significant digit is zero.
  75.     ENDIF
  76.     RETURN
  77.     END
  78. C
  79. C
  80. C
  81. C
  82. C
  83.     FUNCTION ICKCON(IA,N)
  84. C    This routine is used to check if an extended number IA is converging
  85. C    to an integer value. 
  86. C    This routine would probably not be very useful for the average user.
  87. C    The value returned is 1 if IA has converged.  The list significant digit
  88. C    is not checked since it might be just round off.
  89. C    This function return a 1 in IA has converged to an integer otherwise
  90. C    it returns a 0.
  91.     INTEGER*2 IA(*)
  92.     INTEGER*2 N,I
  93.     INTEGER*2 ITEST
  94.     ICKCON=0
  95. C    The number must have already converged past the N-2 th digit before
  96. C    this routine is called.
  97.     ITEST=IA(N-2)
  98.     DO 10 I=3,N-3
  99.       IF(IA(I).NE.ITEST) RETURN
  100. 10    CONTINUE
  101.     ICKCON=1
  102.     RETURN
  103.     END
  104.     SUBROUTINE NEG(IA,N)
  105. C    This routine changes the sign of extended precision number IA.
  106.     INTEGER*2 IA(*),N
  107.     INTEGER*2 I,J
  108.     J=2
  109.     DO 10 I=2,N
  110.       IF(IA(J).NE.0) GOTO 20
  111.       J=J+1
  112. 10    CONTINUE
  113. C    Number was zero
  114.     RETURN
  115. C    The first non-zero digit has been reached
  116. 20    IA(J)=10000-IA(J)
  117.     J=J+1
  118.     DO 30 I=J,N
  119.       IA(I)=9999-IA(I)
  120. 30    CONTINUE
  121.     RETURN
  122.     END
  123. C
  124. C
  125. C
  126. C
  127.     SUBROUTINE PRTEXT(IA,N)
  128. C    This routine writes extended precision numbers to output.
  129. C    This routine could use some work to make a more readable output.
  130.     INTEGER*2 IA(*),N
  131.     INTEGER*2 I,INEG
  132.     INEG=0
  133.     WRITE(*,*) ' '
  134.     IF(IA(N).GT.5000) THEN
  135.       WRITE(*,*) 'Minus'
  136.       CALL NEG(IA,N)
  137.       INEG=1
  138.     ENDIF
  139.     WRITE(*,900) IA(1)
  140. 900    FORMAT(' Exponent  ',I6)
  141.     WRITE(*,910) (IA(I),I=N,2,-1)
  142. 910    FORMAT(10I5)
  143. C    Restore IA to previous value.
  144.     IF(INEG.NE.0) CALL NEG(IA,N)
  145.     RETURN
  146.     END
  147. C
  148. C
  149. C
  150. C
  151.     SUBROUTINE WRTEXT(IUNIT,IA,N)
  152. C    This routine writes extended precision numbers to file IUNIT.
  153. C    This routine could use some work to make a more readable output.
  154.     INTEGER*2 IUNIT,IA(*),N
  155.     INTEGER*2 I,INEG
  156.     INEG=0
  157.     WRITE(IUNIT,*) ' '
  158.     IF(IA(N).GT.5000) THEN
  159.       WRITE(IUNIT,*) 'Minus'
  160.       CALL NEG(IA,N)
  161.       INEG=1
  162.     ENDIF
  163.     WRITE(IUNIT,900) IA(1)
  164. 900    FORMAT(' Exponent  ',I6)
  165.     WRITE(IUNIT,910) (IA(I),I=N,2,-1)
  166. 910    FORMAT(10I5)
  167. C    Restore IA to previous value.
  168.     IF(INEG.NE.0) CALL NEG(IA,N)
  169.     RETURN
  170.     END
  171. C
  172. C
  173. C
  174. C
  175. C
  176.     SUBROUTINE NORMIZ(IA,N)
  177. C    This routine is used by other routines to normalize
  178. C    floating point number after addition and other such processes.
  179.     INTEGER*2 IA(*),N
  180.     INTEGER*2 I,ISHFT
  181.     IF(IA(N).GT.5000) THEN
  182. C      This is a negative number.
  183.       DO 10 I=N,2,-1
  184.         IF(IA(I).NE.9999) GOTO 20
  185. 10      CONTINUE
  186.       I=I+1
  187. 20      ISHFT=N-I
  188.       IF(IA(I).LE.5000) ISHFT=ISHFT-1
  189.     ELSE
  190. C      This is a positive number.
  191.       DO 30 I=N,2,-1
  192.         IF(IA(I).NE.0) GOTO 40
  193. 30      CONTINUE
  194. C      This number is zero
  195.       RETURN
  196. 40      ISHFT=N-I
  197.       IF(IA(I).GE.5000) ISHFT=ISHFT-1
  198.     ENDIF
  199.     IF(ISHFT.EQ.0) RETURN
  200.     IA(1)=IA(1)-ISHFT
  201. C    Shift left the appropriate amount.
  202.     DO 50 I=N,ISHFT+2,-1
  203.       IA(I)=IA(I-ISHFT)
  204. 50    CONTINUE
  205. C    Zero unknown digits.
  206.     DO 60 I=2,ISHFT+1
  207.       IA(I)=0
  208. 60    CONTINUE
  209.     RETURN
  210.     END
  211. C
  212. C
  213. C
  214. C
  215. C
  216.     SUBROUTINE ADD(IA,IB,IC,N)
  217. C    This routine performs extended precision addition.
  218. C    IC = IA + IB. IA and IB are not changed.
  219.     INTEGER*2 IA(*),IB(*),IC(*)
  220.     INTEGER*2 N,IPA,IPB,NML,I
  221.     INTEGER*2 ISGNA,ISGNB,ISGNC,ISGEXT
  222.     INTEGER*2 IT,ICAR
  223.     ISGNA=1
  224.     IF(IA(N).GT.5000) ISGNA=-1
  225.     ISGNB=1
  226.     IF(IB(N).GT.5000) ISGNB=-1
  227. C    Figure out pointers for main loop
  228.     IF(IA(1).GE.IB(1)) THEN
  229.       IC(1)=IA(1)
  230.       IPA=2
  231.       IPB=IA(1)-IB(1)+2
  232.       NML=N-IPB+2
  233.       IF(NML.LT.2) THEN
  234. C        The difference between the two numbers is too large.
  235.         CALL COPY(IA,IC,N)
  236.         RETURN
  237.       ENDIF
  238.       IF(IPB.EQ.2) THEN
  239.         ICAR=0
  240.       ELSE
  241.         ICAR=(IB(IPB-1)+5000)/10000
  242.       ENDIF
  243.     ELSE
  244.       IC(1)=IB(1)
  245.       IPB=2
  246.       IPA=IB(1)-IA(1)+2
  247.       NML=N-IPA+2
  248.       IF(NML.LT.2) THEN
  249. C        The difference between the two numbers is too large.
  250.         CALL COPY(IB,IC,N)
  251.         RETURN
  252.       ENDIF
  253.       ICAR=(IA(IPA-1)+5000)/10000
  254.     ENDIF
  255. C    Main addition loop.
  256.     DO 10 I=2,NML
  257.       IT=IA(IPA)
  258.       IT=IT+IB(IPB)+ICAR
  259.       ICAR=IT/10000
  260.       IC(I)=IT-ICAR*10000
  261.       IPA=IPA+1
  262.       IPB=IPB+1
  263. 10    CONTINUE    
  264. C    Do sign extended part.
  265.     IF(IA(1).GT.IB(1)) THEN
  266.       IF(ISGNB.GE.1) THEN
  267.         ISGEXT=0
  268.       ELSE
  269.         ISGEXT=9999
  270.       ENDIF
  271.       DO 20 I=NML+1,N
  272.         IT=IA(I)
  273.         IT=IT+ISGEXT+ICAR
  274.         ICAR=IT/10000
  275.         IC(I)=IT-ICAR*10000
  276. 20      CONTINUE
  277.     ELSE
  278.       IF(ISGNA.GE.1) THEN
  279.         ISGEXT=0
  280.       ELSE
  281.         ISGEXT=9999
  282.       ENDIF
  283.       DO 30 I=NML+1,N
  284.         IT=IB(I)
  285.         IT=IT+ISGEXT+ICAR
  286.         ICAR=IT/10000
  287.         IC(I)=IT-ICAR*10000
  288. 30      CONTINUE
  289.     ENDIF  
  290. C    Shift right one place if it has overflowed during addition.
  291.     IF(ISGNA*ISGNB.EQ.1) THEN
  292.       ISGNC=1
  293.       IF(IC(N).GE.5000) ISGNC=-1
  294.       IF(ISGNA*ISGNC.EQ.-1) THEN
  295. C        An overflow has occurred.
  296.         IC(1)=IC(1)+1
  297.         ICAR=IC(2)
  298.         DO 40 I=2,N-1
  299.           IC(I)=IC(I+1)
  300. 40        CONTINUE
  301.         IC(N)=ISGEXT
  302. C        Do rounding and return.
  303.         IF(ICAR.GE.5000) CALL INCEX(IC,N)
  304.         RETURN
  305.       ENDIF
  306.     ENDIF
  307. C    Normalize answer.
  308.     CALL NORMIZ(IC,N)
  309.     RETURN
  310.     END
  311. C
  312. C
  313. C
  314. C
  315.     SUBROUTINE MUL(IA,IB,IC,N)
  316. C    This routine multiplies two extended precision numbers.
  317. C    N must be at least 4 or this routine won't work.
  318. C    IC = IA*IB
  319.     INTEGER*2 IA(*),IB(*),IC(*),N
  320.     INTEGER*2 ISNA,ISNB,IXC1,IXC2,ISHFT
  321. C    The computation is carried out with two extra placesIXC1 and IXC2
  322. C    This makes it so precision is not lost if the answer has to be
  323. C    shifted left.  IXC1 is the least significant.
  324.     INTEGER*4 IT,ICAR
  325.     INTEGER*2 I,J,K,IAS
  326.     ISNA=1
  327.     ISNB=1
  328.     IF(IA(N).GT.5000) THEN
  329.       CALL NEG(IA,N)
  330.       ISNA=-1
  331.     ENDIF
  332.     IF(IB(N).GT.5000) THEN
  333.       CALL NEG(IB,N)
  334.       ISNB=-1
  335.     ENDIF
  336.     IC(1)=IA(1)+IB(1)
  337.     IXC1=0
  338.     IXC2=0
  339.     IAS=N+1
  340. C    Step through each IB
  341.     DO 30 I=2,N
  342.       ICAR=0
  343.       IF(I.LE.(N-2)) THEN
  344. C        Round digit that is past least significant.
  345.         IT=IA(IAS-3)
  346.         IT=IT*IB(I)+5000
  347.         ICAR=IT/10000
  348.       ENDIF
  349.       IF(I.LE.(N-1)) THEN
  350. C        Do extended precision 1
  351.         IT=IA(IAS-2)
  352.         IT=IT*IB(I)+ICAR+IXC1
  353.         ICAR=IT/10000
  354.         IXC1=IT-ICAR*10000
  355.       ENDIF
  356. C      Do extended precision 2
  357.       IT=IA(IAS-1)
  358.       IT=IT*IB(I)+ICAR+IXC2
  359.       ICAR=IT/10000
  360.       IXC2=IT-ICAR*10000
  361. C         The first time through the I loop the J loop is not executed
  362.       K=2
  363.       DO 20 J=IAS,N
  364.         IT=IA(J)
  365.         IT=IT*IB(I)+ICAR+IC(K)
  366.         ICAR=IT/10000
  367.         IC(K)=IT-ICAR*10000
  368.         K=K+1
  369. 20      CONTINUE
  370.       IC(K)=ICAR
  371.       IAS=IAS-1
  372. 30    CONTINUE
  373. C    Normalize output if needed rounding when possible.
  374.     IF((IC(N).EQ.0).AND.(IC(N-1).LT.5000)) THEN
  375.       IF((IC(N-1).EQ.0).AND.(IC(N-2).LT.5000)) THEN
  376.         ISHFT=2
  377.       ELSE
  378.         ISHFT=1
  379.       ENDIF
  380.       IC(1)=IC(1)-ISHFT
  381.       DO 40 I=N,ISHFT+2,-1
  382.         IC(I)=IC(I-ISHFT)
  383. 40      CONTINUE
  384.       IF(ISHFT.EQ.2) THEN
  385. C        No rounding is possible in this case.
  386.         IC(3)=IXC2
  387.         IC(2)=IXC1
  388.       ELSE
  389.         IC(2)=IXC2
  390.         IF(IXC1.GE.5000) CALL INCEX(IC,N)
  391.       ENDIF
  392.     ELSE
  393.       IF(IXC2.GE.5000) CALL INCEX(IC,N)
  394.     ENDIF
  395. C    Set sign
  396.     IF(ISNA*ISNB.LT.0) CALL NEG(IC,N)
  397.     RETURN
  398.     END
  399. C
  400. C
  401. C
  402. C
  403.     SUBROUTINE EXSQRT(IA,IB,N,IW)
  404. C    This routine does an extended precision square root.
  405. C    IB = SQRT(IA)
  406. C    It uses an iterative method.  See Arithmetic Operations in Digital
  407. C    Computers by R. K. Richards  D. Van Nostrand Company, Inc. Copyright 1955.
  408. C    B(k+1) = B(k)/2 * (3 - B(k)^2/A)  This series only requires one division.
  409. C    IW is a work area of size 4*N
  410.     INTEGER*2 IA(*),IB(*),IW
  411.     INTEGER*2 N,NUSED,NUSEDL
  412.     INTEGER*2 IALD,IMONE,ITHREE
  413.     DIMENSION IW(N,4)
  414.     REAL*8 R
  415.     IALD=0
  416.     IMONE=-1
  417.     ITHREE=3
  418.     IF(IA(N).GE.5000) THEN
  419.       WRITE(*,*) 'Square root of negative number?'
  420.       STOP
  421.     ENDIF
  422. C    Find first guess.
  423.     CALL CVEXR8(IA,N,R)
  424.     R=SQRT(R)
  425.     NUSED=6
  426.     CALL CVR8EX(R,IB,NUSED)
  427. C    IW(1) is 1/IA at full precision.
  428.     CALL INV(IA,IW(1,1),N,IW(1,2))
  429. C    Beginning of iterative loop
  430. 10    NUSEDL=NUSED
  431.     NUSED=(NUSED-1)*2
  432.     IF(NUSED.GT.N) THEN
  433.       NUSED=N
  434.       IALD=1
  435.     ENDIF
  436.     CALL CHGPRE(IB,NUSEDL,IB,NUSED)
  437.     CALL MUL(IB,IB,IW(1,2),NUSED)
  438. C    Return if convergence is complete
  439.     IF(IALD.EQ.1) THEN
  440.       IF(IEXCMP(IA,IW(1,2),N,IW(1,3)).EQ.0) RETURN
  441.     ENDIF
  442.     CALL CHGPRE(IW(1,1),N,IW(1,3),NUSED)
  443.     CALL MUL(IW(1,2),IW(1,3),IW(1,4),NUSED)
  444.     CALL NEG(IW(1,4),NUSED)
  445.     CALL CVI2EX(ITHREE,IW(1,2),NUSED)
  446.     CALL ADD(IW(1,2),IW(1,4),IW(1,3),NUSED)
  447. C    IW(3) is now equal to   3 - B(k)^2/A
  448.     CALL EXHALF(IW(1,2),NUSED)
  449.     CALL MUL(IW(1,2),IB,IW(1,4),NUSED)
  450.     CALL MUL(IW(1,4),IW(1,3),IB,NUSED)
  451.     GOTO 10
  452.     END
  453. C
  454. C
  455. C
  456. C
  457. C
  458. C
  459.     SUBROUTINE CVI2EX(I2,IA,N)
  460. C    This routine converts an INTEGER*2 number to extended precision.
  461.     INTEGER*2 I2,IA(*)
  462.     INTEGER*2 N,I
  463.     INTEGER*2 ISGN,ITEMP
  464. C    Zero array before start
  465.     DO 10 I=1,N
  466.       IA(I)=0
  467. 10    CONTINUE
  468.     ITEMP=I2
  469.     ISGN=1
  470.     IF(ITEMP.LT.0) THEN
  471.       ISGN=-1
  472.       ITEMP=-ITEMP
  473.     ENDIF
  474.     IF(ITEMP.GE.5000) THEN
  475.       IA(1)=2
  476.       IA(N)=ITEMP/10000
  477.       IA(N-1)=ITEMP-IA(N)*10000
  478.     ELSE
  479.       IA(1)=1
  480.       IA(N)=ITEMP
  481.     ENDIF
  482.       IF(ISGN.EQ.-1) CALL NEG(IA,N)
  483.     RETURN
  484.     END
  485. C
  486. C
  487.     SUBROUTINE COPY(IA,IB,N)
  488. C    This routine copies extended number IA to IB.
  489.     INTEGER*2 IA(*),IB(*),N
  490.     INTEGER*2 I
  491.     DO 10 I=1,N
  492.       IB(I)=IA(I)
  493. 10    CONTINUE
  494.     END
  495. C
  496. C
  497. C
  498. C
  499.     SUBROUTINE EXHALF(IA,N)
  500. C    This routine returns an extended precision 1/2.
  501.     INTEGER*2 IA(*)
  502.     INTEGER*2 N,I
  503.     IA(1)=1
  504.     IA(N)=0
  505.     IA(N-1)=5000
  506.     DO 10 I=2,N-2
  507.       IA(I)=0
  508. 10    CONTINUE
  509.     RETURN
  510.     END
  511. C
  512. C
  513. C
  514. C
  515. C
  516.     SUBROUTINE CVR8EX(R,IA,N)
  517. C    Convert REAL*8 R to an extended precision number IA.
  518. C    N must be greater than or equal to 6.
  519.     REAL*8 R,RT
  520.     INTEGER*2 IA(*)
  521.     INTEGER*2 N,I
  522.     INTEGER*2 ISIGN
  523.     ISIGN=1
  524.     RT=R
  525.     IF(R.LT.0) THEN
  526.       ISIGN=-1
  527.       RT=-RT
  528.     ENDIF
  529.     IA(1)=LOG(RT)/LOG(10000.)
  530.     RT=RT/(10000.**IA(1))
  531.     IA(1)=IA(1)+1
  532.     IF(RT.GE.5000) THEN
  533.       IA(1)=IA(1)+1
  534.       RT=RT/10000.
  535.     ENDIF
  536.     DO 10 I=N,N-4,-1
  537.       IA(I)=RT
  538.       RT=RT-IA(I)
  539.       RT=RT*10000
  540. 10    CONTINUE
  541.     DO 20 I=N-5,2,-1
  542.       IA(I)=0
  543. 20    CONTINUE
  544.     IF(ISIGN.LT.0) CALL NEG(IA,N)
  545.     RETURN
  546.     END
  547. C
  548. C
  549. C
  550.     SUBROUTINE CVEXR8(IA,N,R)
  551. C    Convert an extended precision number IA to REAL*8 R.
  552. C    N must be greater than or equal to 6.
  553.     REAL*8 R,RT,EXP
  554.     INTEGER*2 IA(*)
  555.     INTEGER*2 N,I
  556.     INTEGER*2 ISIGN
  557.     ISIGN=1
  558.     IF(IA(N).GE.5000) THEN
  559.       ISIGN=-1
  560.       CALL NEG(IA,N)
  561.     ENDIF
  562.     RT=0
  563.     EXP=1
  564.     EXP=(EXP/10000)**5
  565.     DO 10 I=N-4,N
  566.       RT=RT+IA(I)*EXP
  567.       EXP=EXP*10000
  568. 10    CONTINUE
  569.     EXP=10000
  570.     R=RT*(EXP**IA(1))
  571.     IF(ISIGN.LE.0) R=-R
  572.     RETURN
  573.     END
  574. C
  575. C
  576. C
  577. C
  578.     SUBROUTINE CHGPRE(IA,NA,IB,NB)
  579. C    This routine changes the precision of an extended precision number.
  580. C    This routine is useful for many iterative solutions where the full
  581. C    precision is not needed until the final iteration.
  582. C    IA and IB can be the same array.
  583. C    Input is IA of size NA, output is IB of size NB.
  584.     INTEGER*2 IA(*),IB(*),ICAR
  585.     INTEGER*2 NA,NB,I,J
  586.     IB(1)=IA(1)
  587.     IF(NA.GT.NB) THEN
  588. C      IA has the most precision.
  589.       ICAR=0
  590.       J=NA-NB+1
  591.       IF(IA(J).GT.5000) ICAR=1
  592.       DO 10 I=2,NB
  593.         J=J+1
  594.         IB(I)=IA(J)
  595. 10      CONTINUE
  596.       IF(ICAR.EQ.1) CALL INCEX(IB,NB)
  597.     ELSE
  598. C      IB has the most precision.
  599.       J=NB
  600.       DO 20 I=NA,2,-1
  601.         IB(J)=IA(I)
  602.         J=J-1
  603. 20      CONTINUE
  604.       DO 30 I=2,J
  605.         IB(I)=0
  606. 30      CONTINUE
  607.     ENDIF
  608.     RETURN
  609.     END
  610. C
  611. C
  612. C
  613. C
  614. C
  615.     SUBROUTINE INV(IX,IY,N,IW)
  616. C    This routine does an extended precision inverse.
  617. C    IY=1/IX
  618. C    It uses an iterative method.  A routine more like MUL would probably 
  619. C    be more efficient and would certainly take less storage, but I haven't
  620. C    gotten around to writing it yet.
  621. C    Y(I+1)=Y(I)*(2-X*Y(I))
  622. C    IW is a work area of size 3*N
  623.     INTEGER*2 IX(*),IY(*),IW
  624.     INTEGER*2 N,NUSED,NUSEDL
  625.     INTEGER*2 IALD,ISGN,INEGT
  626.     DIMENSION IW(N,3)
  627.     REAL*8 R
  628.     INEGT=-2
  629.     IALD=0
  630.     ISGN=1
  631.     IF(IX(N).GE.5000) THEN
  632.       CALL NEG(IX,N)
  633.       ISGN=-1
  634.     ENDIF
  635. C    Find first guess.
  636.     CALL CVEXR8(IX,N,R)
  637.     R=1/R
  638.     NUSED=6
  639.     CALL CVR8EX(R,IY,NUSED)
  640. 10    NUSEDL=NUSED
  641.     NUSED=(NUSED-1)*2
  642.     IF(NUSED.GT.N) THEN
  643.       NUSED=N
  644.       IALD=1
  645.     ENDIF
  646.     CALL CHGPRE(IX,N,IW(1,1),NUSED)
  647.     CALL CHGPRE(IY,NUSEDL,IY,NUSED)
  648.     CALL MUL(IW(1,1),IY,IW(1,2),NUSED)
  649. C    Return if convergence is complete
  650.     IF(IALD.EQ.1) THEN
  651.       IF(ICKCON(IW(1,2),N).EQ.1) THEN
  652.         IF(ISGN.LT.0) CALL NEG(IY,N)
  653.         RETURN
  654.       ENDIF
  655.     ENDIF
  656.     CALL CVI2EX(INEGT,IW(1,1),NUSED)
  657.     CALL ADD(IW(1,1),IW(1,2),IW(1,3),NUSED)
  658.     CALL NEG(IW(1,3),NUSED)
  659.     CALL MUL(IY,IW(1,3),IW(1,1),NUSED)
  660.     CALL COPY(IW(1,1),IY,NUSED)
  661.     GOTO 10
  662.     END
  663. C
  664. C
  665. C
  666. C
  667. C
  668. C
  669. C
  670.     FUNCTION IEXCMP(IA,IB,N,IW)
  671. C    This routine compares two extended precision numbers.
  672. C    The returned value is 1 if IA > IB, 0 if IA is within roundoff of IB, etc.
  673. C    IW is a work area of size N.
  674. C    Only works for positive numbers
  675.     INTEGER*2 IA(*),IB(*),IW(*)
  676.     INTEGER*2 N,I
  677.     REAL*8 A,B
  678.     REAL DIF
  679. C    First check if the numbers are close.
  680.     CALL CVEXR8(IA,N,A)
  681.     CALL CVEXR8(IB,N,B)
  682.     DIF=(A-B)/(A+B)
  683.     IF(ABS(DIF).GT.1E-12) THEN
  684.       IEXCMP=1
  685.       IF(DIF.LT.0) IEXCMP=-1
  686.       RETURN
  687.     ENDIF
  688. C    The values are close
  689.     CALL NEG(IB,N)
  690.     CALL ADD(IA,IB,IW,N)
  691.     CALL NEG(IB,N)
  692. C    Check for zero
  693.     IF(IA(1)-IW(1).GE.N-2) THEN
  694. C      Within roundoff of zero
  695.       IEXCMP=0
  696.       RETURN
  697.     ENDIF
  698.     DO 10 I=2,N
  699.       IF(IW(I).NE.0) GOTO 20
  700. 10    CONTINUE
  701.     IEXCMP=0
  702.     RETURN
  703. C    Results are not zero
  704. 20    IEXCMP=1
  705.     IF(IW(N).GT.5000) IEXCMP=-1
  706.     RETURN
  707.     END
  708. 
  709.